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Abstract — In this paper we propose a model-order reduction 
method for chemical reaction networks governed by general 
enzyme kinetics, including the mass-action and Michaelis- 
Menten kinetics. The model-order reduction method is based 
on the Kron reduction of the weighted Laplacian matrix which 
describes the graph structure of complexes in the chemical 
reaction network. We apply our method to a yeast glycolysis 
model, where the simulation result shows that the transient 
behaviour of a number of key metabolites of the reduced-order 
model is in good agreement with those of the full-order model. 



I. Motivation and introduction 

The dynamics of a biochemical reaction network can 
be described by a set of differential equations involving 
a stoichiometric matrix and reaction rates. The solution to 
these equations describes the evolution and dynamics of the 
concentrations of all the metabolites or biochemical species 
in the network. The stoichiometric matrix, which consists 
of (positive and negative) integer elements, captures the 
basic conservation laws of the reactions. The reaction rates 
are functions of the concentrations and they prescribe the 
dynamics that take place at each reaction in the network 
based on the underlying enzyme kinetics. The common type 
of reaction rates are mass-action kinetics and Michaelis- 
Menten-type kinetics laws which can include competitive and 
non-competitive inhibition mechanism. 

Thus a biochemical reaction network £ can be written as 

/ x = Sv(x) + S b v b 

: \ y =Cx 

where x and v b denote the vectors of metabolite concentra- 
tions and boundary fluxes, respectively, y is the measured 
concentrations (or the metabolites of interest), the function 
v is the reaction rates and C denotes the observation matrix 
of appropriate dimension. The matrix S and S b are the 
stoichiometric matrices of the internal and boundary fluxes, 
respectively. 

In this paper, we describe a new model-order reduction 
method to £ by reducing the number of reactions and 
metabolites that are involved, such that the dynamics of 
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the measured concentration of the reduced model is close 
to the full-order one (with the external fluxes). Our model- 
order reduction method is based on the Kron reduction of 
the underlying weighted Laplacian describing graph structure 
of complexes in the chemical reaction networks. A similar 
technique is also employed in the Kron reduction method for 
model reduction of resistive electrical networks described in 
[6]. 

In our previous works [9], [13], we have proposed a 
framework for describing the dynamics of balanced chemical 
reaction networks governed by mass-action kinetics where, 
using the notion of complexes, the dynamics can be writ- 
ten by a complex stoichiometric matrix and a symmetric 
weighted Laplacian matrix. The symmetric Laplacian matrix 
corresponds to the graph structure of the reactions and the 
weights come from the underlying equilibrium constants 
in each reaction. The model-order reduction method that 
we proposed in [9], [13] is based on the application of 
Kron reduction to this Laplacian matrix. In this paper, we 
extend the method to a general class of biochemical reaction 
networks. 

For biochemical reaction networks, the model-order reduc- 
tion technique for nonlinear systems is still underdeveloped. 
The singular perturbation method and quasi steady-state 
approximation (QSSA) approach are the most commonly 
used techniques, where the reduced state contains a part of 
the metabolites of the full model. In the thesis by Hardin [3], 
the QSSA approach is extended by considering the higher- 
order approximation in the computation of quasi steady- 
state. Sunnaker et al. in [11] proposed a reduction method 
by identifying variables that can be lumped together and 
can be used to infer back the original state. In Prescott 
& Papachristodoulou [8], a method to compute the upper- 
bound of the error is proposed based on sum-of-squares 
programming. The application of these techniques to general 
kinetics laws, such as Michaelis-Menten, poses a significant 
computational problem. 

Our proposed method, on the other hand, offers a simple 
approach to model-order reduction that can be extended 
directly to general kinetics. Moreover, the resulting reduced- 
order model retains the structure of the kinetics and gives 
result to a reduced biochemical reaction network, which 
enables a direct biochemical interpretation. 

We show the application of our model reduction technique 
to the yeast glycolysis model described in [14]. We have 
simulated the transient behaviour of the metabolites that are 
not eliminated during the model reduction procedure. We 
show that there is a good agreement between the transient 
behaviour of the concentration of most of such metabolites 



when comparing the full network to the reduced network. 

The paper is organized as follows. In Section II, we intro- 
duce tools from stoichiometry of reactions and graph theory 
that are required to derive our mathematical formulation. In 
section III, we explain enzyme kinetics, and then derive our 
formulation. In section IV, we propose our model reduction 
method. In section V, we show the application of our method 
to a yeast glycolysis model and in section VI, we present 
conclusions based on our results. 

Notation: The space of n dimensional real vectors is denoted 
by R n , and the space of m x n real matrices by R mxn . 
The space of n dimensional real vectors consisting of all 
strictly positive entries is denoted by R™ and the space of n 
dimensional real vectors consisting of all nonnegative entries 
is denoted by Wf , The rank of a real matrix A is denoted 
by rank A. t m denotes a vector of dimension m with all 
entries equal to 1. The time-derivative 4f (£) of a vector x 
depending on time t will be usually denoted by x. I n denotes 
an identity matrix of dimension n. 

Define the mapping Ln : R™ — > R m , x t-¥ Ln(a;), as 
the mapping whose z-th component is given as (Ln(.T)) i :— 
ln(xi). Similarly, define the mapping Exp : R m — > 
R™, x H> Exp(ir), as the mapping whose 7-th component 
is given as (Exp(x)) i := exp(a;j). 

II. Chemical reaction network structure 

In this section, we introduce the tools necessary in order 
to derive our mathematical formulation of the dynamics of 
chemical reaction networks. First we introduce the concept 
of stoichiometric matrix of a reaction network. We then 
introduce the concept of a complex graph, which was first 
introduced in the work of Horn & Jackson and Feinberg ([4], 
[5], [2]). 

A. Stoichiometry 

Consider a chemical reaction network involving m chemi- 
cal species (metabolites), among which r chemical reactions 
take place. The basic structure underlying the dynamics of 
the vector x S R+ of concentrations Xi,i — 1, • • • , to, of 
the chemical species is given by the balance laws x = Sv, 
where S is an to x r matrix, called the stoichiometric matrix. 
The elements of the vector v € R r are commonly called the 
(reaction) fluxes or rates. The stoichiometric matrix S, which 
consists of (positive and negative) integer elements, captures 
the basic conservation laws of the reactions. It already 
contains useful information about the network dynamics, 
independent of the precise form of the reaction rate v(x). 
Note that the reaction rate depends on the governing law 
prescribing the dynamics of the reaction network. 

We now show how to construct the stoichiometric matrix 
for a reaction network with the help of an example shown 
in Fig. [1] 

X 4 



Note that the above reaction has 5 reactions among four 
species X\, X 2 , X3 and X4. Since the stoichiometric matrix 
maps the space of reactions to the space of species, it has 
dimension 4x5. The entry of S corresponding to the 7 th 
row and j th column is obtained by subtracting the number 
of moles of the 7 th species on the product side from that on 
the substrate side for the j th reaction. Thus 
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for the reaction network depicted in Fig. [T] 

B. Complex graph 

The structure of a chemical reaction network cannot be 
directly captured by an ordinary graph. Instead, we will 
follow an approach originating in the work of Horn and 
Jackson [4], introducing the space of complexes. The set of 
complexes of a chemical reaction network is simply defined 
as the union of all the different left- and righthand sides 
(substrates and products) of the reactions in the network. 
Thus, the complexes corresponding to the network ([T]i are 
X 1 + 2X 2 , X 3 , 2X 1 + X 2 and X 4 . 

The expression of the complexes in terms of the chemical 
species is formalized by an to x c matrix Z, whose a-th 
column captures the expression of the a-th complex in the 
to chemical species. For example, for the network depicted 
in Figure [T] 

10 2 
2 10 
10 
1 
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Xi + 2X 2 ^ X 3 — 2Xi 

Fig. 1. Example of a reaction network 



X<y 



We will call Z the complex stoichiometric matrix of the 
network. Note that by definition all elements of the matrix 
Z are non-negative integers. 

Since the complexes are the left- and righthand sides of the 
reactions, they can be naturally associated with the vertices of 
a directed graph Q with edges corresponding to the reactions. 
Formally, the reaction a — > /3 between the a-th (reactant) 
and the /3-th (product) complex defines a directed edge with 
tail vertex being the a-th complex and head vertex being the 
/3-th complex. The resulting graph will be called the complex 
graph. 

Recall, see e.g. [1], that any graph is defined by its 
incidence matrix B. This is a c x r matrix, c being the 
number of vertices and r being the number of edges, with 
(a,j)-th element equal to —1 if vertex a is the tail vertex 
of edge j and 1 if vertex a is the head vertex of edge j, 
while otherwise. Thus the structure of the complex graph 
is described by a c x r incidence matrix B. 

Obviously, there is a close relation between the matrix Z 
and the stoichiometric matrix S. In fact, it is easily checked 
that 

S = ZB, hence x = ZBv(x) 
with B the incidence matrix of the complex graph. 



In many cases of interest, especially in biochemical reac- 
tion networks, chemical reaction networks are intrinsically 
open, in the sense that there is a continuous exchange 
with the environment (in particular, inflow and outflow of 
chemical species and connection to other reaction networks). 
In this paper, we are particularly interested in inflows to and 
outflows from the complexes of the network. This will be 
modeled by a vector w& <E E c consisting of c boundary (or, 
exchange) fluxes, leading to an extended model 



x 



ZBv(x) + Zvi, 



(1) 



III. The dynamics of networks governed by 

ENZYME KINETICS 

Recall that for a chemical reaction network, the relation 
between the reaction rates and species concentrations de- 
pends on the governing laws of the reactions involved in the 
network. In this section, we explain this relation for reaction 
networks governed by enzyme kinetics. In other words, if 
v denotes the vector of reaction rates and x denotes the 
species concentration vector, we show how to construct v(x) 
for reaction networks governed by enzyme kinetics. 

Let Zs denote the column of the complex stoichiometric 
matrix Z corresponding to the substrate Sj of the j-th 
reaction of a chemical reaction network. Let kj denote the 
rate constant of the f h reaction of the network. Then the 
reaction rate of the j th reaction of the chemical reaction 
network between the j th substrate Sj and the j th product 
Vj is given by 



Vj(x) — dj (x)kj exp (zJ j .Ln(x)), 



(2) 



where for j = 1, . . . , r, dj : — > R+ is a rational function 
of its argument. Note that if the governing law of the j th 
reaction is mass action kinetics, then dj(x) = 1. 
As an example, consider the reaction 



Xx + 3X 2 



X 3 + 3X 4 



(3) 



governed by Michaelis Menten kinetics which is the most 
common form of enzyme kinetics. For i = 1, . . . , 4, let a;, 
denote the concentration of the species X,. Define x := 
\x\ X2 X3 X4] T . The matrices B, Z and S for the 
reaction ([3]l are given by 



B = 



Z = 



"1 


0" 
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s = 


-3 
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1 
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3 



Note that Z s is given by Z s = [l 3 0] Let K x , 
K 2 , -K3 and denote the "Michaelis" constants of the 
species Xi, X 2 , X 3 and ^^respectively. Let Vt denote the 
maximum reaction rate of ( 3 1. In this case, k 
net rate of the reaction ([3]) is given by 



v f The 



where the expression for d(x) depends on the model used to 
define the rate of the reaction. One possibility for d(x) is 

1 



d(x) = 



1 1 

1 + Kt 



3x 2 
K 2 



K 1 



X3_ 



3X4 



d[x)kx\x\ = d(x)kexp (Z];Ln(x)) 



Based on the formulation in Q, we can describe the 
complete reaction network dynamics as follows. Let the 
reaction rate for the complete set of reactions be given 
by the vector v(x) = \vi(x) ■■■ v r (x)\ . For every 
(7,71 £ {1, ••• , c}, define a„(i) := kjdj{x) if (c, 7r) = 
(Sj,Vj), j £ {1, . . . , r} and a a7r :— otherwise. Define the 
weighted adjacency matrix A of the complex graph as the 
matrix with (er, 7r)-th element a CT7r , where a, w £ {1, • • • , c}. 
Furthermore, define L(x) := A(.t) — A{x), where A is the 
diagonal matrix whose (p, p)-th element is equal to the sum 
of the elements of the p-th column of A. Hereafter we refer 
to L(x) as the weighted Laplacian of the complex graph 
associated with the given network with species concentration 
vector x. It is a matter of straightforward verification to check 
that t^L(x) = 0. It can also be verified that the vector 
Bv(x) for the mass action reaction rate vector v(x) is equal 
to — L(x)Exp (Z T Ln(a;)), where the mapping Exp : M c — > 
M.S_ has been defined at the end of the Introduction. Hence 
the dynamics can be compactly written as 

x = -ZL(x)Exp (Z T Ln(x)) 

A similar expression of the dynamics corresponding to mass 
action kinetics, in less explicit form, was already obtained 
in [10]. 

A. The linkage classes of a complex graph 

A linkage class of a chemical reaction network is a maxi- 
mal set of complexes {C\, . . . ,Ck} such that Ci is connected 
by reactions to Cj for every i,j £ {1, ...,&}, i 7^ j. It 
can be easily verified that the number of linkage classes (£) 
of a network, which is equal to the number of connected 
components of the complex graph corresponding to the 
network, is given by I = c— rank(£?) (the number of linkage 
classes in the terminology of [4], [2]). 

IV. Model reduction 

For many purposes one may wish to reduce the number of 
dynamical equations of a chemical reaction network in such 
a way that the behaviour of a number of key metabolites 
is approximated in a satisfactory way. We propose a novel 
method for model reduction of chemical reaction networks 
governed by enzyme kinetics. Our method is based on the 
Kron reduction method for model reduction of resistive 
electrical networks described in [6]; see also [12]. 

A. Description of the method 

Our model reduction method is based on reduction of the 
complex graph associated with the network. It is based on the 
following result regarding Schur complements of weighted 
Laplacian matrices. 

Proposition 4.1: Consider a chemical reaction network 
with weighted Laplacian matrix L(x) £ R cxc corresponding 



L{x) 



(4) 



to the concentration vector x. Let V denote the set of vertices 
of the complex graph associated with the network. Then 
for any subset of vertices V r C V, the Schur complement 
L(x) of L(x) with respect to the indices corresponding to 
V r satisfies the following properties: 

1) All diagonal elements of L(x) are positive and off- 
diagonal elements are nonnegative for x € K™. 

2) tjL(x) = 0, where c := c - dim(V r ). 

Proof: (7.) Follows from the proof of [7, Theorem 
3.11]; see also [12] for the case of symmetric L. 

(2.) Without loss of generality, assume that the last c — c 
rows and columns of L(x) correspond to the vertex set V r . 
Consider a partition of L(x) given by 

Ln(x) L 12 {x) 
L 2 i(x) L 22 (x) 

where L lx (x) G R £x£ , L 12 (x) G R £x ( c - £ \ L 21 {x) G 
K ( c -c)xc and i 22 (a;) e R (c-e)x(c-e)_ It is easy t0 see that 

L(ar) = Lii(x) - Li 2 (a;)i22(a;) _:L i2i(a;) 
Since t^L(x) — 0, we obtain 

tjLnix) + t^_ £ L 21 {x) = tjL 12 {x) + l^_ a L 22 (x) = 
Using the above equations, we get 

lTL(x) = lJ(L u (x) - L 12 {x)L 22 (x)- l L 2l (x)) 
= -\ T c _- c L 21 {x) + \ T c ^ c L 22 {x)L 22 (x)- x L 21 (x) = 

■ 

From the above result, it follows that L(x) obeys all the 
properties of the weighted Laplacian of a reaction network 
corresponding to the complex graph with vertex set V — V r . 
Thus Proposition |4.1| can be directly applied to the complex 
graph, yielding a reduction of the chemical reaction network 
by reducing the number of complexes. Consider a reaction 
network with boundary fluxes described as in equation ([TJ 

S: x = Z (v b - L{x)Exp(Z T Ln(x))) (5) 

Reduction will be performed by deleting certain complexes 
in the complex graph, resulting in a reduced complex graph 
with weighted Laplacian L(x). Furthermore, leaving out 
the corresponding columns of the complex stoichiometric 
matrix Z one obtains a reduced complex stoichiometric 
matrix Z (with as many columns as the remaining num- 
ber (c) of complexes in the complex graph). Consider a 
corresponding partition of L given by equation Q. Define 
P := \lg —L X2 L 22 ] ■ Now consider the reduced model 



E : x = Z [Pv b - L(x)Exp(Z T Ln(x)) 

= ZP(v h - L{x)Exp(Z T Ln(x))) . (6) 

Note that S is again a chemical reaction network governed 
by enzyme kinetics, with a reduced number of complexes 
and with, in general, a different set of boundary fluxes and 
reactions (edges of the complex graph). 

An interpretation of the reduced network E can be given 
as follows. Consider a subset V r of the set of all complexes, 



which we wish to leave out in the reduced network. Consider 
the partition of L[x) as given by equation Q and corre- 
sponding partitions of Z and v b given by 



z = [z x z 2 ] 



Vb 



[vbi 



Vb2\ 



where V r corresponds to the last part of the indices (denoted 
by 2), in order to write out the dynamics of E as 



Vbl 


- Z 


v b2 _ 





Lu(x) L 12 (x) 
L 2 i(x) L 22 (x) 



Exp(ZfLn(») 
Exp(ZjLnO)) 



Consider now the auxiliary dynamical system 



m 




Vbl 




Ln(x) 


Lx 2 (x) 




Wi 


2/2 




Vb2 




L 2 i(x) 


L 2 2(x)_ 




W 2 



where we impose the constraint y 2 = 0. It follows that 

w 2 = -L 2 2(xy 1 (v b 2 - L 2 i(x)wi), 
leading to the reduced auxiliary dynamics 

Vi = (vbi - Li2{x)L 2 2(xy 1 v b 2) - L(x)wi 
= Pvb — L(x)wi 

Putting back in wi = Exp(Zj r Ln(x)), making use of x — 
Z\i)\ + Z2V2 = Ziyx_= Zyi, we then obtain the reduced 
network E given in (|6jl. An appropriate choice of V r will 
ensure that some of the elements of x have derivative zero 
in E leading to lesser number of state variables in E as 
compared to E. 

Assume that a given biochemical reaction network is 
asymptotically stable around an equilibrium. When perturbed 
from the equilibrium, such a reaction network has certain 
species reaching their equilibrium much faster than the 
remaining ones. The principle behind our model reduction 
method is to impose the condition that complexes entirely 
made up of such species remain at constant concentrations. 

Example 4.2: We now consider an example of a simple 
reversible reaction network governed by Michaelis-Menten 
kinetics and apply the model reduction procedure described 
in this paper. This reaction is shown below. 



X\ + X2 ^- X3 + X4 ^- X5 



Xa 



(7) 



For i — 1, ...,6, let Xi denote the concentration of Xi. 
Let v\ and V2 denote overall reaction rates in the forward 
direction of the first and the second reversible reactions 
respectively. For i = 1, ... ,4, let K mi i denote the Michaelis 
constant of Xi for the first reaction. For i — 3, . . . , 6, let 
K m2 i denote the Michaelis constant of Xi for the second 
reaction. For i = 1,2, let k{ and k\ denote the forward 
and reverse rate constants of the i th reaction. Define x := 
col(xi,x 2 ,x 3 ,X4,,x 5 ,xe), 



il mil lv rni1 

p 2 (x) := 1 + — + 



1 



X 3 Xi 



K m i3 -"-mi 4 



Then v x = hj^l. 



and V2 = 



consider the reduced network 



Now 



X* 



obtained by deleting the complex X% + X4 from the network 
(|7j. Applying the procedure described in this section, we 
obtain the following expression for the overall rate (v) in the 
forward direction of the reduced network: 



1 



K„ 



K„ 



+ 



(8) 



K„ 



where k[, k%, K m3 x, K m3 2, K m35 , K m3 e are positive con- 
stants. Note that equation <|8j represents the rate of a reaction 
governed by Michaelis-Menten kinetics. Thus for this exam- 
ple, our procedure yields a reduced model with 6 parameters 
while the original model has 12 parameters. 

B. Effect of Model Reduction 

In this section, we show the effect of our model reduction 
method on a particular type of networks with a single linkage 
class. In other words, we give an interpretation of our 
reduced model in terms of its corresponding full model for 
a particular type of networks. Note that deletion of a set of 
complexes in one linkage class does not affect the reactions 
of the other linkage classes of the network. 



Full Network: C x ^ C 2 ^ C 3 
Reduced Network: C\ ^ C 3 



C 



(9) 
(10) 



Consider a reaction network with reversible reactions 
occuring between consecutive elements of the set of distinct 
complexes {C\,C2, ■ ■ ■ ,Cn} as in d9V The reduced network 



obtained by deleting the complex C2 looks as in ( 10 1. Instead 
of the two reversible reactions, C\ ^ C2 and C2 C3 in the 
full network, there is one reaction C\ ^ C3 in the reduced 
network. This reaction is a reversible reaction governed by 
enzyme kinetics. The expression for the rate of this reaction 
in terms of the expression for the rates of the first two 
reversible reactions of the full network can be found using 
the technique described in this section. All the remaining 
reactions of the reduced network occur in the same way as 
in the full network. 

Assuming that the network is asymptotically stable around 
a certain equilibrium, when perturbed from the equilibrium, 
the transient behaviour of the metabolites involved in the 
complexes of the reduced model will approximately be the 
same as that of the full model if the metabolites involved in 
C2 reach their steady states much faster than the remaining 
metabolites. In this case, we can safely impose the condition 
that the metabolites involved in C 2 are at constant concentra- 
tion in order to obtain the reduced model ( fT0| ) with similar 
transient behaviour as that of (|9j. The rule of induction may 
be applied in order to further reduce the model by deleting 
more complexes. 

A special case of networks d9| is the following: 




Fig. 2. The yeast glycolysis model as discussed in [14]. The figure is taken 
from [14]. 



Deletion of the complex C2 in this case is equivalent to 
deletion of the linkage class from the network. Such a dele- 
tion provides a close approximation to the original network 



(11) 



if the reaction (Hi has very little effect on the dynamics 
of the network, i.e if the reaction reaches its steady state 
much faster than the remaining reactions of the network, 
assuming the network is asymptotically stable around a 
certain equilibrium. 

Observe that for the model reduction of a given asymp- 
totically stable biochemical reaction network governed by 
enzyme kinetics, it is important to determine which of the 
complexes are to be deleted so that the reduced model 
approximates the full model reasonably well. 

V. Yeast Glycolysis Model 

We have applied our model reduction procedure on the 
model of yeast glycolysis described in [14]. The schematic 
of the model is shown in Fig. [2] The corresponding detailed 
mathematical model can be found in [14]. 

For the modelling, we have considered the following as 
external fluxes: 

• in/efflux of extracellular glucose (Glco) to/from intra- 
cellular glucose (Glci); 

• constant influx of trehalose to Glci; 

• constant efflux of trehalose from glucose 6-phosphate 
(G6P); 

• constant efflux of glycerol from dihydroxyacetone phos- 
phate (DHAP); 
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Fig. 4. Concentration of Glci vs time 
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Fig. 3. The reduced model of yeast glycolysis. The figure is a modified 
form of the original glycolysis model in [14]. 
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• efflux of succinate from pyruvate (PYR) proportional to 
the concentration of PYR; 

• efflux of acetate from acetaldehyde (AcAld) propor- 
tional to the concentration of AcAld; 

• in/efflux of ethanol to/from AcAld catalyzed by alcohol 
dehydrogenase (ADH). 

The reactions of the network are governed by enzyme 
kinetics and we can write the equations of the model in 
the same form as |5]). Conditions of a glucose pulse as 
described in [14] are assumed for the model. According to 
these conditions, for t < 0, concentrations of Glco and ATP 
are 0.2 mM and 5 mM respectively, and for t > 0, these are 
equal to 5 mM and 2.5 mM respectively. It is assumed that 
ethanol (EtoH) has a constant concentration and the network 
is at equilibrium for t < 0. With these conditions, it is found 
that the model is asymptotically stable as t — > oo. 

Various combinations of complexes have been consid- 
ered for deletion in order to obtain a reduced model 
that closely mimics the transient behaviour of the original 
model. It is found that deletion of the complexes G6P, 
3-phosphoglycerate (3PG), 2-phosphoglycerate (2PG) and 
phosphoenoylpyruvate (PEP) from the original model pro- 
duces the best results in this respect. The schematic of the 
reduced model is shown in Fig. [3] 



i 1 1 1 1 1 
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Fig. 5. Concentration of PYR vs time 

It is observed that there is a good agreement between the 
transient behaviours of most of the metabolites when com- 
paring the original model to the reduced model. This implies 
that in the original model, under the given conditions of the 
glucose pulse, the metabolites G6P, 3PG, 2PG and PEP reach 
their steady states much faster than the remaining metabolites 
of the network. Thus imposing that these metabolites stay at 
a constant concentration has very little effect on the dynamics 
of the network. The graphs of comparison of the evolution 
of concentrations of Glci and PYR are shown in Figures [4] 
and [5] respectively. 

VI. Conclusion 

In this paper, we have outlined a method for model reduc- 
tion of biochemical reaction networks that are asymptotically 
stable around an equilibrium. When perturbed from the equi- 
librium, such a reaction network has certain species reaching 
their equilibrium much faster than the remaining ones. The 
principle behind our model reduction method is to impose 



the condition that a few of such species remain at constant 
concentrations. This is achieved by deletion of complexes 
entirely made up of such species from the complex graph 
of the network. Apart from a reduction in the number of 
state variables, our model reduction method also leads to a 
reduction in the number of reactions and parameters of the 
model. We have applied our method on a yeast glycolysis 
model and have observed a good agreement between the 
transient behaviours of the original and the reduced models. 
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